c                                                                               
c *** program filter -- applies the cubic smoothing spline to a time series.                            
c
      integer ymax
      parameter (nmax=2500,ymax=6000,iu14=14,iu15=15)                                                                               
      dimension tr(nmax),sp(nmax),res(nmax),cl(nmax),wrk(nmax),wrk1(nmax),wrk2(nmax)
      real*8 arps(nmax),phi(nmax),pm,wk1(nmax),wk2(nmax),wkm(nmax)
c
      integer iyr(ymax),lyr(ymax),nyr(ymax),ncol
      integer p,i,ii
      real*8 val(nmax,ymax)
      character isame,norm
      character*20 tits(ymax)
c
      character fil1*80,fil2*80,frm*40,tit*80      
      logical frun
c
c February 7,2008 -pjk
c New in version ARS41d is the addition of a call to PLPLOT_fonts2.c to set an environment 
c vairable to include path to /Applications/Absoft/EXAMPLES/PLPlot/examples the folder that 
c contains those six pesty font files for plplot graphics
c
      call PLPLOT_fonts( )    !new for Absoft10.1 OS X 10.5 Intel
c
      write(*,'(/1x,''Enter the input file (/ to exit)'',t45,''==> '',$)')
      read(*,'(a)')fil1
      if(fil1(1:1).eq.'/')goto 200
      open (iu14,file=fil1,access='sequential',status='old')             
      write(*,'(1x,''Enter the output file'',t45,''==> '',$)')
      read(*,'(a)')fil2
      open (iu15,file=fil2,access='sequential',status='unknown')             
      call rdtabdat(iu14,nmax,ymax,val,iyr,lyr,nyr,ncol,tits,sentin)
      close(iu14)
c
      write(*,'(/1x,''Enter the 50% frequency response cutoff'',t45,''==> '',$)')
      read(*,*)ls 
      norm=' '
      write(*,'(1x,''Normalize the input data first? y/n'',t45,''==> '',$)')
      read(*,'(a)')norm      
      if(norm.eq.'Y')norm='y'
      if(norm.eq.'N')norm='n'
      isame=' '
      write(*,'(1x,''Maximum x-axis range for all plots? y/n'',t45,''==> '',$)')
      read(*,'(a)')isame
      if(isame.eq.'Y')isame='y'
      if(isame.eq.'N')isame='n'

    5 write(*,'(1x,''What column to start plots? Range 1 -'',i4,t45,''==> '',$)')ncol
      read(*,*)istcol
      if(istcol.lt.1.or.istcol.gt.ncol)goto 5
    6 write(*,'(1x,''What column to end plots?   Range 1 -'',i4,t45,''==> '',$)')ncol
      read(*,*)iencol
      if(istcol.lt.1.or.istcol.gt.ncol)goto 6
      if(iencol.le.istcol)goto 6
c
ccc      do p=1,ncol  !start a big loop here to process each series in the file
      do p=istcol,ncol  !start a big loop here to process each series in the file
c
        ifyr=iyr(p);ilyr=lyr(p);n=nyr(p)
        nsave=n
        ifyrsave=ifyr
        ilyrsave=ilyr
  
        nskp1=0
        nskp2=0
        do i=1,n
          datum=val(i,p)
          if(datum.eq.sentin)then
            nskp1=nskp1+1
          elseif(datum.ne.sentin)then
            goto 10
          endif
        enddo
   10   do i=1,n
          ii=n-i+1
          datum=val(ii,p)
          if(datum.eq.sentin)then
            nskp2=nskp2+1
          elseif(datum.ne.sentin)then
            goto 15
          endif
        enddo
   15   n=n-nskp1-nskp2
        iyr(p)=iyr(p)+nskp1
        ifyr=iyr(p)
        lyr(p)=lyr(p)-nskp2
        nyr(p)=n
        do i=1,n
          ii=i+nskp1
          cl(i)=val(ii,p)
          tr(i)=cl(i)
        enddo
        np=0
        do i=1,n
          datum=cl(i)
          if(datum.ne.sentin)then
            np=np+1
            wrk1(np)=cl(i)
          endif
        enddo
        call mnsd(wrk1,np,xm,sd)
        do i=1,n
          datum=tr(i)
          if(datum.eq.sentin)then
            if(norm.eq.'y')then
              tr(i)=0.0
              cl(i)=0.0
            elseif(norm.eq.'n')then
              tr(i)=xm
              cl(i)=xm
            endif
          else
            if(norm.eq.'y')then
              tr(i)=(tr(i)-xm)/sd
              cl(i)=0.0
            elseif(norm.eq.'n')then
              tr(i)=tr(i)-xm
              cl(i)=xm
            endif
          endif
        enddo
        if(norm.eq.'y')xm=0.0
        if(norm.eq.'n')xm=0.0

c
c *** calculate prediction error filter coefficients
c
ccc        ip=ls
ccc        lg=300
ccc        iopt=0
ccc        write(*,'(/1x,''*** Computing prediction error filter'',
ccc     *    '' (length '',i3,'') for extrapolation''/)')ip
ccc        call memcof(tr,n,ip,pm,arps,wk1,wk2,wkm)
ccc        write(*,'(1x,''ip: '',i3)')ip
ccc        write(*,'(1x,10f8.3)')(arps(j),j=1,ip)
c
c *** extend the series backwards and forwards using pef estimates from mempr
c
ccc        next=ls
ccc        np=n+2*next
ccc        do i=1,nmax
ccc          wrk(i)=0.0
ccc        enddo
ccc        do i=1,n
ccc          ii=i+next
ccc          wrk(ii)=tr(i)
ccc        enddo
ccc        call extend(wrk,n,next,arps,ip,0)
ccc        call spline(np,wrk,sp,pp,ls)         
        call spline(n,tr,sp,pp,ls)         
ccc        do i=1,n
ccc          ii=i+next
ccc          tr(i)=wrk(ii)+xm
ccc          sp(i)=sp(ii)+xm
ccc        enddo
        write(iu15,'(/1x,''Input series is: '',a)')tits(p)
        write(iu15,'('' year  unfilt  spline  resids'')')
        do i=1,n
          res(i)=wrk(i)-sp(i)
          write(iu15,'(i5,3f8.3)')i+ifyr-1,tr(i),sp(i),res(i)
        enddo
        if(isame.eq.'n')then        
          call plot1(nmax,tr,sp,cl,nyr(p),iyr(p),tits(p),norm)
        elseif(isame.eq.'y')then
          do i=1,nmax
            wrk1(i)=0.0
            wrk2(i)=0.0
            cl(i)=xm
          enddo
          do i=1,n
            ii=i+nskp1
            wrk1(ii)=tr(i)
            wrk2(ii)=sp(i)
          enddo
          if(nskp1.gt.0)then
            do i=1,nskp1
              wrk1(i)=xm
              wrk2(i)=xm
            enddo 
          endif
          call plot1(nmax,wrk1,wrk2,cl,nsave,ifyrsave,tits(p),norm)
        endif
      enddo
      close(iu15)
  200 stop                                                                      
      end                                                                       
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine rdtabdat(iu,nmax,ymax,val,iyr,lyr,nyr,ncol,tits,sentin)
c
	integer nmax,ymax,iu,ncol
	integer iyr(ymax),lyr(ymax),nyr(ymax)
	real*8 val(nmax,ymax)
	character*20 tits(ymax)
c
	integer p,j,k,nrow
	real   sentin
c
	Write(*,'(1x,''Enter number of columns including years'',t45,''==> '',$)')
	read(*,*)ncol
	write(*,'(1x,''Enter missing value flag, e.g., -999.00'',t45,''==> '',$)')
	read(*,*)sentin
c
	read(iu,*)(tits(p),p=1,ncol)
	nrow=1
10	read(iu,*,end=20)(val(nrow,k),k=1,ncol)
	   nrow=nrow+1
	goto 10
20	nrow=nrow-1
c
	do k=2,ncol
	  iyr(k-1)=int(val(1,1))    ! by defulat set individual series' first year to the youngest posible date in file
	  lyr(k-1)=int(val(nrow,1)) ! by defulat set individual series' last year to the oldest posible date in file
ccc	  do p=1,nrow  !test to see if iyr and lyr are other than first or last values
	  do p=2,nrow-1  !test to see if iyr and lyr are other than first or last values
	    if( (val(p,k).ne.sentin) .and. (val(p-1,k).eq.sentin) )iyr(k-1)=int(val(p,1))
	    if( (val(p,k).ne.sentin) .and. (val(p+1,k).eq.sentin) )lyr(k-1)=int(val(p,1))
	  enddo
	  nyr(k-1)=lyr(k-1)-iyr(k-1)+1
	enddo
c
	ncol=ncol-1
c
	do k=1,ncol
	  j=1
	  do p=1,nrow
	    if(val(p,k+1).ne.sentin)then
	       val(j,k)=val(p,k+1)
	       j=j+1
	    elseif( (val(p,k+1).eq.sentin) .and. (p.gt.nyr(k)) )then
	       val(p,k)=0.0
	    
	    endif
	  enddo
	enddo
c
	do p=1,ncol
	    tits(p)=tits(p+1)
	enddo
	
c
	return
	end
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
	subroutine plot1(nmax,wrk,sp,cl,ny,iy,tit,norm)
c
c	  write(winu9,'(/13x,''COLOR TABLE ______________________________________________'')')
c	  write(winu9,'(13x,''White(0) Red(1)  Yellow(2)  Green(3)  Aquamarine(4) Pink(5) '')')
c	  write(winu9,'(13x,''Wheat(6)  Grey(7)  Browm(8)  Blue(9)  BlueViolet(10)'')')
c	  write(winu9,'(13x,''Cyan(11)  Turquoise(12)  Magenta(13)  Salmon(14)  Black(15)'')')
c	  write(winu9,'(13x,''___________________________________________________________'',/)')

c
c	argument
	  integer nmax,ny,iy
	  dimension sp(nmax),wrk(nmax),cl(nmax)
	  character norm
	  character*20 tit
c	local
	   integer p,red,green,blue
	   real*4 x(10000),mny,mxy,mx,mnx	
	   logical bg
c
	mny=100000.0;mxy=-100000.0
	red=255;green=255;blue=255
	bg=.false.
c
	do p=1,ny
	  if(wrk(p).gt.mxy)mxy=wrk(p)
	  if(sp(p).gt.mxy)mxy=sp(p)
	  if(wrk(p).lt.mny)mny=wrk(p)
	  if(sp(p).lt.mny)mny=sp(p)
	  x(p)=iy+p-1
	enddo
c	
	mnx=iy;mx=iy+ny-1
	write(*,*)mxy,mny
	write(*,*)mx,mnx

ccc       mxy=6.0
ccc       mny=-6.0

c	
  	call plsdev('Mac1')
c
c set background color
	if(.not.bg)then
	  call plscol0(0,255,255,255)              !change default color map black>white 
	  call plscol0(15,0,0,0)                   !change default color map white>black
	  call plscolbg(red,green,blue)	        !change background color if other than white is desired
	elseif(bg)then
	  call plscol0(15,255,255,255)              !change default color map white>black
	  call plscol0(0,0,0,0)                   !change default color map black>white
	  call plscolbg(red,green,blue)	        !change background color if other than white is desired
	endif  
c
c 	initialize plplot
	call plinit()
c 	set lable color
	call plcol(15)   !color box and lables
c 	define the viewport one column by one rows
	call plssub(1,1)
	call plenv(mnx,mx,mny,mxy,0,2)
	if(norm.eq.'y')then
  	  call pllab('Years','Standard Normal Deviates',tit)
	else
  	  call pllab('Years','Original Data Units',tit)
	endif
c 	set line style
	call pllsty(1)
c 	set curve width
	call plwid(1)
c 	set line color
	call plcol(1)  !color lines
	call plline(ny,x,wrk)            !plot wrk
c 	set curve width
	call plwid(2)
c 	set curve color
	call plcol(9)  !color lines
	call plline(ny,x,sp)	          !plot spline curve
c 	set curve width
	call plwid(2)
c	set lable color
	call plcol(15)   !color box and lables
	call plline(ny,x,cl)	          !plot mean curve
c	
	call plend()
c	
	return
	end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                           
      subroutine rank(data,nobs,nmax)                                              
      dimension data(nmax)                                                         
      r=nobs                                                                    
      do 10 i=1,nobs                                                            
      big=data(i)                                                               
      save=data(i)                                                              
      k=i                                                                       
      do 20 j=i,nobs                                                            
      if(data(j).le.big) go to 20                                               
      big=data(j)                                                               
      k=j                                                                       
   20 continue                                                                  
      data(k)=save                                                              
      data(i)=big                                                               
   10 continue                                                                  
      return                                                                    
      end
c  
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      subroutine extend(y,n,ncast,arps,ip,iopt)
      real*4 y(1)
      real*8 arps(1),pf,pb,zero
      data zero/0.0d0/
ccc      write(2,'(10f8.4)')(arps(j),j=1,ip)			!debug
      do 15 i=1,ncast
        pf=zero
        pb=zero
        nf=ncast+n+i
        nb=ncast+1-i
        do 10 j=1,ip
          if(iopt.eq.0.or.iopt.eq.+1)pf=pf+arps(j)*dble(y(nf-j))
          if(iopt.eq.0.or.iopt.eq.-1)pb=pb+arps(j)*dble(y(nb+j))
   10   continue
        y(nf)=sngl(pf)
        y(nb)=sngl(pb)
ccc        write(2,'(1x,2i5,f9.5,i4,f9.5)')i,nb,y(nb),nf,y(nf)	!debug
   15 continue
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
*sbr spline: fits cubic spline curve to time series
*
*  derived from imsl library routines by edward r cook,
*  lamont-doherty geological observatory, and modified by
*  richard l holmes, tree-ring laboratory, university of arizona.
*
*  modified by r l holmes 10 jun 1988
*  modified by e r cook to full double precision 4 apr 1991
*
********************************************************************
*
* calculates a cubic spline function to filter a time series.
* 'stiffness' of spline is specified by parameter 'ls'
*
* parameters:
*
* n:  number of values in time series
* y:  time series array to be modeled
* f:  array of cubic spline function values computed
* p:  value computed in this routine, used to calculate spline
* ls: stiffness (50% frequency response) of spline
* a:  working array of dimension (n,4)
*
********************************************************************
*
      subroutine spline(n,y,f,pp,ls)
      parameter(mxy=10000)
      implicit real*8 (a-h,o-z)
      dimension a(mxy,4),c1(4),c2(3)
      real*4 y(n),f(n),pp
      data c1 /1.d0,-4.d0,6.d0,-2.d0/
      data c2 /0.d0,.33333333333333d0,1.33333333333333d0/
      data pct,pi / .50d0, 3.1415926535897935d0 /
      if(n .lt. 4)then
      f(1)=-9998.
      return
      endif
      if(n.gt.mxy)then
        write(*,'(/t2,''*** series too long for spline subroutine ***''/
     *           t2,'' *** enlarge mxy parameter in subroutine ***''/)')
        stop
      endif
      nm2=n-2
      v=dble(ls)
      arg=(2.d0*pi)/v
      p=(6.d0*(dcos(arg)-1.d0)**2.d0)/(dcos(arg)+2.d0)
      pp=sngl(p)
      do 1 i=1,nm2
      do 1 j=1,3
      a(i,j)=c1(j)+p*c2(j)
    1 a(i,4)=dble(y(i))+c1(4)*dble(y(i+1))+dble(y(i+2))
      a(1,1)=c2(1)
      a(1,2)=c2(1)
      a(2,1)=c2(1)
      nc=2
* begin ludapb
      rn=dble(1.d0/(dfloat(nm2)*16.d0))
      d1=1.d0
      d2=0.d0
      ncp1=nc+1
      if(nc .ne. 0)then
* initialize zero elements
      do 5 i=1,nc
      do 5 j=i,nc
      k=ncp1-j
    5 a(i,k)=0.d0
* 'i' is row index of element being computed
* 'j' is column index of element being computed
* 'l' is row index of previously computed vector
*       being used to compute inner product
* 'm' is column index
      endif
      do 60 i=1,nm2
      imncp1=i-ncp1
      i1=max0(1,1-imncp1)
      do 60 j=i1,ncp1
      l=imncp1+j
      i2=ncp1-j
      sum=a(i,j)
      jm1=j-1
      if(jm1 .gt. 0)then
      do 25 k=1,jm1
      m=i2+k
   25 sum=sum-a(i,k)*a(l,m)
      endif
      if(j .eq. ncp1)then
      if(a(i,j)+sum*rn .le. a(i,j))then
      f(1)=-9999.
      write(*,'(''>> sbr spline: matrix not positive definite'')')
      return
      endif
      a(i,j)=1.d0/dsqrt(sum)
* update determinant
      d1=d1*sum
   35 if(dabs(d1) .gt. 1.d0)then
      d1=d1*.0625d0
      d2=d2+4.d0
      goto 35
      endif
   47 if(dabs(d1) .le. .0625d0)then
      d1=d1*16.d0
      d2=d2-4.d0
      goto 47
      else
      goto 60
      endif
      endif
      a(i,j)=sum*a(l,ncp1)
   60 continue
* end ludapb / begin luelpb
* solution ly = b
      nc1=nc+1
      iw=0
      l=0
      do 15 i=1,nm2
      sum=a(i,4)
      if(nc .gt. 0)then
      if(iw .ne. 0)then
      l=l+1
      if(l .gt. nc) l=nc
      k=nc1-l
      kl=i-l
      do 45 j=k,nc
      sum=sum-a(kl,4)*a(i,j)
   45 kl=kl+1
      else
      if(sum .ne. 0.d0)iw=1
      endif
      endif
   15 a(i,4)=sum*a(i,nc1)
* solution ux = y
      a(nm2,4)=a(nm2,4)*a(nm2,nc1)
      n1=nm2+1
      do 30 i=2,nm2
      k=n1-i
      sum=a(k,4)
      if(nc .gt. 0)then
      kl=k+1
      k1=min0(nm2,k+nc)
      l=1
      do 50 j=kl,k1
      sum=sum-a(j,4)*a(j,nc1-l)
   50 l=l+1
      endif
   30 a(k,4)=sum*a(k,nc1)
* end luelpb
      do 2 i=3,nm2
    2 f(i)=a(i-2,4)+c1(4)*a(i-1,4)+a(i,4)
      f(1)=a(1,4)
      f(2)=c1(4)*a(1,4)+a(2,4)
      f(n-1)=a(nm2-1,4)+c1(4)*a(nm2,4)
      f(n)=a(nm2,4)
      do 3 i=1,n
    3 f(i)=y(i)-f(i)
      return
      end
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
csbr memcof
c
      subroutine memcof(data,n,m,pm,cof,wk1,wk2,wkm)
      implicit real*8 (a-h,o-z)
      real*4 data(1)
      real*8 cof(1),wk1(1),wk2(1),wkm(1)
      data zero/0.0d0/,one/1.0d0/,two/2.0d0/
      rn=n
      p=zero
      do 11 j=1,n
        p=p+dble(data(j))*dble(data(j))
11    continue
      pm=p/rn
      wk1(1)=dble(data(1))
      wk2(n-1)=dble(data(n))
      do 12 j=2,n-1
        wk1(j)=dble(data(j))
        wk2(j-1)=dble(data(j))
12    continue
      do 17 k=1,m
        pneum=zero
        denom=zero
        do 13 j=1,n-k
          pneum=pneum+wk1(j)*wk2(j)
          denom=denom+wk1(j)**2+wk2(j)**2
13      continue
        cof(k)=two*pneum/denom
        pm=pm*(one-cof(k)**2)
        if(k.ne.1)then
          do 14 i=1,k-1
            cof(i)=wkm(i)-cof(k)*wkm(k-i)
14        continue
        endif
        if(k.eq.m)return
        do 15 i=1,k
          wkm(i)=cof(i)
15      continue
        do 16 j=1,n-k-1
          wk1(j)=wk1(j)-wkm(k)*wk2(j)
          wk2(j)=wk2(j+1)-wkm(k)*wk1(j+1)
16      continue
17    continue
      pause 'never get here'
      end
*sbr mnsd
*
      subroutine mnsd(x,n,xm,sd)
      real*4 x(n),xm,sd,sq,sm,rn
      xm=0.0
      sd=0.0
c
c compute the mean and variance of the series
c
      sq=0.0
      sm=0.0
      do 10 i=1,n
        sm=sm+x(i)
   10 continue
      rn=n
      xm=sm/rn
      do 15 i=1,n
        sq=sq+(x(i)-xm)**2
   15 continue
      sd=sqrt(sq/(rn-1.0))
      return
      end